#include <iostream>
#include <vector>
#include <string>
#include <fstream>
#include <map>
#include <stdlib.h>
#include <lineages.h>
#include <algorithm>
#include <stdio.h>
#include <sstream>
#include <math.h>
#include <time.h>
#include <shared.h>

using namespace std;

int fitnessType;
int robustType;
int disturbType;

double baseDelRate = 1.0;

void usage()
{
	cerr << "Generate a set of random lineages that is based on the a real lineage" << endl;
	cerr << "This program is specifically designed for H.roretzi lineage," << endl;
	cerr << "where the two half-lineages are exactly the same, so that the random lineages retain this feature." << endl;
	cerr << "Usage: mirrorRandom func lin num seed collect" << endl;
	cerr << "Example: mirrorRandom 1 hro.alm 100 123 hro.random.lineages" << endl;
	cerr << "Note: All the lineages other than the H.roretzi lineage is handled differently, see \"mkRandom\" instead" << endl;
	cerr << "Parameters (all are needed) :"<< endl;
	cerr << " func      : The method of generating the random lineages" << endl;
	cerr << "             1. total random coalescence" << endl;
	cerr << "             2. keep the shape of the lineage tree, but randomly shuffle the terminal cells" << endl;
	cerr << "             3. generate \"num\" simulated trees for each of the maximum depths in M to (2M + 2 - X) ," << endl;
	cerr << "                where X is the maximum depth of the real lineage," << endl;
	cerr << "                and M is the theoretically lowest possible maximum depth of the real lineage." << endl;
	cerr << "             4. generate random lineages so that the collection of terminal depths in each of them" << endl;
	cerr << "                is exactly the same as the real lineage" << endl;
	cerr << "             5. generate random lineages so that the collection of terminal depths in each of them" << endl;
	cerr << "                is a bootstrap sample of the collection of terminal depths in the real lineage" << endl;
	cerr << "             6. shuffle the terminal cells within their own depth, do not change lineage tree shape" << endl;
	cerr << " lin       : Name of a lineage file (one cell per line)" << endl;
	cerr << "             The random lineages are generated based on this lineage." << endl;
	cerr << "             The file name must contain the string \"hro\" to indicate it is a H.roretzi lineage." << endl;
	cerr << "             The program need to know what lineage is being used because" << endl;
	cerr << "             it need to \"understand\" the names of functional classification" << endl;
	cerr << "             I recommend only those lineage files downloadable from the program's google code project site." << endl;
	cerr << " num       : Generate how many simulated lineages." << endl;
	cerr << "             Note that for \"func\" = 3, there will be a multiplication of \"num\" lineages." << endl;
	cerr << " seed      : Seed for randomization" << endl;
	cerr << " collect   : File name for the lineage collection (one lineage per line)" << endl;
	cerr << "             It must follow the same rule as \"lin\"" << endl;
	cerr << " " << endl;
	return;
}

int main(int argc, char *argv[])
{
    if(argc != 6)
	{
		usage();
		exit(255);
	}
    int func = atoi(argv[1]);
    string realLineagesFile = argv[2];
    int simLoop = atoi(argv[3]);
    int randSeed = atoi(argv[4]);
    string outputFile = argv[5];
    srand(randSeed);

    vector<int> termTypes;
    vector<string> termLineages;

    loadRealLineages(realLineagesFile,termTypes,termLineages);
    Lineages realLin(termTypes,termLineages);
    realLin.getDepthStat();

    int outputBuffer = 100;
    vector<vector<int> > typesPool;
    vector<vector<string> > lineagesPool; // buffer for output

    vector<int> halfTypes;
    vector<string> halfLineages;
    demirror(termLineages,termTypes,halfLineages,halfTypes); // starting from these half lineages

	vector<string> randHalfLineages; // resulting half lineage
    vector<int> randTypes; // resulting full lineage
	vector<string> randLineages; // resulting full lineage

    if(func == 1) // random topology
    {
        for(int i=0; i<simLoop; i++)
        {
            generateRandomLineageTopology(halfTypes.size(),randHalfLineages);
            mirror(randHalfLineages,halfTypes,randLineages,randTypes);
            typesPool.push_back(randTypes);
            lineagesPool.push_back(randLineages);
            if(i == 0)
            {
                writeLineageCollection(outputFile,typesPool,lineagesPool);
                typesPool.clear();
                lineagesPool.clear();
            }
            else if((i % outputBuffer == 0) or (i == simLoop - 1) )
            {
                appendLineageCollection(outputFile,typesPool,lineagesPool);
                typesPool.clear();
                lineagesPool.clear();
            }
        }
    }
    else if(func == 2) // random order of terminal cells
    {
        for(int i=0; i<simLoop; i++)
        {
            vector<int> tmpTypes = halfTypes;
            random_shuffle(tmpTypes.begin(),tmpTypes.end());
            mirror(halfLineages,tmpTypes,randLineages,randTypes);
            typesPool.push_back(randTypes);
            lineagesPool.push_back(randLineages);

            if(i == 0)
            {
                writeLineageCollection(outputFile,typesPool,lineagesPool);
                typesPool.clear();
                lineagesPool.clear();
            }
            else if((i % outputBuffer == 0) or (i == simLoop - 1) )
            {
                appendLineageCollection(outputFile,typesPool,lineagesPool);
                typesPool.clear();
                lineagesPool.clear();
            }
        }
    }
    else if(func == 3) // random topology with controlled depth limit
    {
        int minDepth,maxDepth;
        minDepth = (int)(log(termTypes.size()) / log(2)) + 1;
        maxDepth = (abs(minDepth - (int)realLin.depthMax) + 1) * 2 + (int)realLin.depthMax;
        for(int j=minDepth; j <= maxDepth; j++)
        {
            for(int i=0; i<simLoop; i++)
            {
                generateRandomLineageTopologyLimitDepth(halfTypes.size(),j-1,randHalfLineages);
                // half lineage, so depth - 1
				mirror(randHalfLineages,halfTypes,randLineages,randTypes);
                typesPool.push_back(randTypes);
                lineagesPool.push_back(randLineages);
                if((i == 0) and (j == minDepth))
                {
                    writeLineageCollection(outputFile,typesPool,lineagesPool);
                    typesPool.clear();
                    lineagesPool.clear();
                }
                else if((i % outputBuffer == 0)  or (i == simLoop - 1) )
                {
                    appendLineageCollection(outputFile,typesPool,lineagesPool);
                    typesPool.clear();
                    lineagesPool.clear();
                }
            }
        }
    }
    else if(func == 4) // random topology with same depths of all terminal cells
    {
        vector<int> halfDepths;
        for(int i=0; i<simLoop; i++)
        {
            halfDepths.clear();
            for(unsigned int j=0;j<halfLineages.size();j++)
            {
                halfDepths.push_back(halfLineages[j].length());
            }
            random_shuffle(halfDepths.begin(),halfDepths.end());
            depthFixedCoalescent2(halfDepths,randHalfLineages);
			mirror(randHalfLineages,halfTypes,randLineages,randTypes);
            typesPool.push_back(randTypes);
            lineagesPool.push_back(randLineages);
            if(i == 0)
            {
                writeLineageCollection(outputFile,typesPool,lineagesPool);
                typesPool.clear();
                lineagesPool.clear();
            }
            else if( (i % outputBuffer == 0) or (i == simLoop - 1) )
            {
                appendLineageCollection(outputFile,typesPool,lineagesPool);
                typesPool.clear();
                lineagesPool.clear();
            }
        }
    }
    else if(func == 5) // randomly sample real depths then build tree for them
    {
    	vector<int> sampledDepths;
    	vector<string> randLineages;
		for(int i=0; i<simLoop; i++)
        {
        	sampledDepths.clear();
			do{
			    sampledDepths.clear();
				for(unsigned int j=0;j<halfLineages.size();j++)
				{
					sampledDepths.push_back(halfLineages[(int)random(0.0,(double)halfLineages.size())].length());
				}
			} while(!goodDepthSet(sampledDepths));

            depthFixedCoalescent2(sampledDepths,randHalfLineages);
			mirror(randHalfLineages,halfTypes,randLineages,randTypes);
            typesPool.push_back(randTypes);
            lineagesPool.push_back(randLineages);
            if(i == 0)
            {
                writeLineageCollection(outputFile,typesPool,lineagesPool);
                typesPool.clear();
                lineagesPool.clear();
            }
            else if( (i % outputBuffer == 0) or (i == simLoop - 1) )
            {
                appendLineageCollection(outputFile,typesPool,lineagesPool);
                typesPool.clear();
                lineagesPool.clear();
            }
        }

    }
    else if(func == 6) // shuffled termini within their depth
    {
        for(int i=0; i<simLoop; i++)
        {
            depthFixedShuffle(halfLineages,randHalfLineages);
            mirror(randHalfLineages,halfTypes,randLineages,randTypes);
            typesPool.push_back(randTypes);
            lineagesPool.push_back(randLineages);

            if(i == 0)
            {
                writeLineageCollection(outputFile,typesPool,lineagesPool);
                typesPool.clear();
                lineagesPool.clear();
            }
            else if((i % outputBuffer == 0) or (i == simLoop - 1) )
            {
                appendLineageCollection(outputFile,typesPool,lineagesPool);
                typesPool.clear();
                lineagesPool.clear();
            }
        }
    }

    return 0;
}
